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ABSTRACT 

A radiant heat flux gage calibration system exists in the 
Flight Loads Laboratory at NASA's Dryden Flight Research Center. 
This calibration system must be well understood if the heat flux gages 
calibrated in it are to provide useful data during radiant heating 
ground tests or flight tests of high speed aerospace vehicles. A part of 
the calibration system characterization process is to develop a 
numerical model of the flat plate heater element and heat flux gage, 
which will help identify errors due to convection, heater element 
erosion, and other factors. 

A 2-dimensional mathematical model of the gage-plate 
system has been developed to simulate the combined problem 
involving convection, radiation and mass loss by chemical reaction. A 
fourth order finite difference scheme is used to solve the steady state 
governing equations and determine the temperature distribution in the 
gage and plate, incident heat flux on the gage face, and flat plate 
erosion. Initial gage heat flux predictions from the model are found to 
be within ±7% of experimental results, 

NOMENCLATURE 

A = area 

Bh = blowing parameter 

Cp = specific heat 

D = diameter of the heat flux gage 

DI = channel diameter 

D i2 = oxygen diffusion coefficient 

F = shape factor 

Gr = Grashof number 

g = acceleration due to gravity 

H = height 

h t = average natural convection heat transfer coefficient 

I = electrical current through the graphite plate 

k = thermal conductivity 

L = length 

Q = heat generation per unit volume of graphite plate 

q = heat flux 

m, c =5 oxygen mass concentration near the graphite plate surface 

m 1(j = oxygen mass concentration in the ambient air 

m = mass transfer rate from horizontal graphite plate surface 

n u = oxygen mass flux near the graphite plate surface 


n 2s = carbon dioxide mass flux near the graphite plate surface 
Pr = Prandtl number 

R = electrical resistance 

Re = Reynold number 

Sc = Schmidt number 

T = temperature 

AT = temperature difference between the surface and the 

ambient 

TH = thickness 
Greek letters: 

P lra = mass transfer driving force 

E = radiation emissivity 

i] = mass transfer conductance 

rj t . = mass transfer conductance 

X = specific resistance 

v = kinematic viscosity 

p = density of the graphite plate 

a = Stefan-Boltzmann constant 

tp = Prandtl number function 

Subscripts: 

1 = oxygen 

2 = carbon dioxide 

a = ambient air 

g = thermal gage 

p = graphite plate 

s = surface 

x = horizontal direction 

y as vertical direction 

water = water cooled 

Superscripts: 

0 = zero mass transfer 

d = heated plate facing down 

u = heated plate facing up 



1. INTRODUCTION 

Uncertainties in radiant heat flux gage calibrations are typically 
assumed to be on the order of ±10%. The National Institute of 
Standards and Technology (NIST) is currently working to 
significantly reduce this uncertainty at heat flux levels up to 200 
kW/m 2 as described in HoImberg(1997); Grosshandler(1997); 
Murthy(1997). However, aerodynamic heating in hypersonic flight 
can generate heat flux levels well in excess of this level. The Flight 
Loads Laboratory (FLL) at NASA’s Dryden Flight Research Center 
performs radiant thermal-structural tests on aerospace vehicle 
structures which represent aerodynamic heat fluxes up to 1,100 
kW/m 2 . This paper describes a numerical analysis which is part of an 
effort undertaken in the FLL to reduce radiant heat flux calibration 
uncertainties at heat fluxes up to 1 ,100 kW/m 2 , 

A radiant heat flux gage calibration system in the FLL uses an 
electrically heated graphite plate as a heat source (A detailed 
description is given by hom(1993)). The main objective of the project 
described in this paper is to numerically model the flat plate heater 
and compare numerical and experimental results in an effort to 
quantify and, if possible, reduce errors in the calibration system. An 
initial 2-dimensional model of the system has been developed and the 
results have been compared to data from experiment. The numerical 
model provides the temperature distribution in both the graphite plate 
and the heat flux gage as well as the heat flux at the gage face. The 
finite difference method was used in the numerical model. The 
boundary conditions on the front face of the heat flux gage include 
heat received by radiation from the graphite plate and lost to the 
ambient atmosphere by natural convection. Heat transfer by forced 
convection of water is modeled on the back face of the gage. Heat is 
transferred from the graphite plate by natural convection and radiation 
to the ambient surroundings as well as the heat flux gage on one side 
of the plate. Also, oxidation of the graphite plate, which impacts 
radiation view factors, is included in the model. 

2. EXPERIMENT 



Figure 1: Heat Flux Calibration System, Side View 


The heat flux gage calibration system in the Flight Loads 
Laboratory at NASA’s Dryden Flight Research Center consists of a 
graphite plate clamped between two copper electrodes (see Fig. 1). 
The graphite plate undergoes ohmic healing as low voltage (<24 
volts), high current rup to 2400 Amps), electric power is passed 
through the plate. The plate can reach temperatures up to 2800° C. 


Calibrations can be transferred from one heat flux gage to another by 
placing one heat flux gage on each side of the graphite plate, as shown 
in Fig. 2. The heat generation is assumed to be uniform. 



Figure 2: Heat Flux Calibration System, Top View 


The experiments performed to support the characterization of the 



Data Acquisition 
System 


Figure 3: Characterization Test System Schematic 

flat plate heater utilized a circular foil heat flux gage to measure heat 
flux. The calibration supplied by the gage manufacturer was used. 
The second heat flux gage was eliminated. Other measurements 
include surface temperature, voltage across the plate, current through 
the plate, and plate erosion. 

An infra-red pyrometer, described by Matthews(1989), measured 
plate surface temperature at the center of the clear side of the plate 
(Fig. 3). The target spot size was approximately 1 cm diameter. The 
pyrometer uses a laser based technique to measure emissivity and 
compute a temperature corrected for emissivity. It also incorporates 
through-the-lens sighting to aid in aiming and focusing the pyrometer 
optics on the measurement location. A wire mesh grid placed over 
the flat graphite plate provided features on which to aim and focus. 
The grid was removed prior to testing. 

Voltage was measured at the copper electrodes. Electric current 
passing through the plate was measured as well. The voltage, current, 
and heat flux data were acquired on a personal computer based data 
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acquisition system using voltmeters, a GPIB interface, and associated 
data acquisition software. Temperature data was logged every 10 
seconds on board the infra-red pyrometer and printed out after each 
test. Plate thickness was measured at the center of the plate using a 
micrometer before and after each test. Erosion of the top and bottom 
of the plates were measured post-test using a taper gage and flat 
surface. 

Tests were performed at nominal heat fluxes of 100 kW/m 2 , 200 
kW/m\ 350 kW/m\ and 450 kW/m 2 . A test was run for 2, 4, 6, 8, and 
10 minutes for each nominal heat flux. This matrix allowed the 
acquisition of voltage, current, and plate erosion data as a function of 
heat flux (i.e., plate temperature) and run time The experimental data 
presented in this paper were obtained during the nominal 450 kW/m 2 , 
10 minute run. 

A typical test proceeds as follows: 

1 . Measure thickness of a new graphite plate 

2. Install the plate between the copper electrodes. The four set 
screws visible in Fig. 1 secure the plate and are uniformly torqued 
using a torque limiting screwdriver. 

3. Apply minimum power to the flat plate and visually check 
for uniform heating as the plate begins to glow. 

4. Remove power from plate. Retorquc set screws and 
recheck plate heating or proceed to next step as required. 

5. Start data acquisition and recording 

6. Apply power to the plate by setting the manual power 
supply control to the predetermined level which corresponds to the 
desired nominal heat flux. 

7. Continue power application for the desired length of time 
without readjusting the manual power supply control. 

8. Remove power from the plate when desired run time has 
been reached. 

9. Remove graphite plate from between the electrodes. 

10. Measure eroded plate thickness as well as top and bottom 
erosion. 

Some additional remarks are necessary regarding steps 3, 6 and 7 
above. Some small amount of erosion occurs when the plate is 
checked for uniform heating (step 3 above). The minimum power 
possible is used for this run and the total energized lime is less than i 
minute. The erosion caused in step 3 is extremely small and cannot 
be measured without disturbing the test setup, which would introduce 
still more errors into the heat flux measurement. 

The nominal heat flux values mentioned earlier are only target 
values. The resolution of the manual power supply control allows the 
target heat flux to be set within 30 kW/m 2 in step 6. For this reason, 
the heat flux data presented in this report will differ slightly from the 
nominal, or target, value. 

Step 7, above, indicates that the manual power supply control is 
not adjusted once it is set for the desired heat flux. The constant 
power supply setting and plate erosion result in significant decrease in 
the measured heat flux (up to 20%) and flat plate current (up to 32%) 
while causing a voltage increase of up to 5%. 

3. NUMERICAL MODELING 

The numerical model developed in this study includes 
conduction heat transfer within the graphite plate and heat flux gage 
model, heat generation within the graphite plate, and convective, 
radiative, and mass transfer boundary conditions. A schematic 
diagram of the experiment set-up is shown in Fig. 4. The numerical 
model represents a vertical plane which passes through the flat plate 
and heat flux gage midway between the electrodes (the hatched lines). 
Figure 5. shows the configuration used in this study to simulate the 2- 
Dimensional model. This two dimensional model is being used to 


develop modeling techniques which incorporate the various physical 
phenomena and will be extended to three dimensions in the near 
future. 

The problem is assumed to be quasi-steady state. That is, heat 
conduction within the graphite plate and heat flux gage is essentially 
steady state at any given moment in time even though the boundary 
conditions (i.e., current through the plate, mass loss, etc.) vary slowly 
with time. This assumption is implemented in the numerical model 
by using an implicit finite difference method to solve the two 
dimensional, steady state heat conduction equations at discrete time 
steps (1 second intervals). The mass loss from the graphite plate over 
1 second is computed after each time step. The next time step 
conduction solution proceeds with new plate thickness and height, as 
determined from the previous mass loss calculation, and new current 
through the flat plate, as measured during the experiment. 



Figure 4: Schmematic of Flat Plate and Heat Flux Gage With 
The Computation Cross Section 
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Figure 5: Computation Regime 


3.1 Numerical model for the graphite plate 

As shown in Fig. 5, the length of the plate is L p , the height is H r , 
and the thickness is TH p . The steady state governing energy equation 
for the graphite plate is: 
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where k p is the graphite plate thermal conductivity. Information 

provided by the graphite manufacturer indicates that the type ATJ 
graphite used in the plate is relatively isotropic. Therefore, the value 
used for k p is the same in both directions. The thermal conductivity 

of graphite steadily decreases from room temperature at 25 C to 
3500° C. The functional dependence is provided by the graphite plate 
vendor. The temperature dependence used in the numerical 
calculation is interpolated from the curve supplied by the vendor. The 
graphite plate is heated by passing electric current through it. The 
model incorporates this as a heat generation term, Q m , as 


i 2 -r 

L p H p TH p 


(W hi f) 


(2) 


where I is the electrical current and R is the electrical resistance in 
the direction of current flow. The resistance R is calculated as 


R = 


H p TH p 


(3) 


x = TH p ; 

-k p ^- = q p+ h cy (T-T e )-m y AH c (7) 

At y=0, and y=H p heat is lost by radiation and natural convection 
and is generated by combustion. Eq. (4) and (5) represent this balance. 
At x=0, heat is lost by the radiation and convection to the ambient 
surroundings and gained because of the chemical reaction. At x=TH p , 
the boundary conditions include radiation and convection to the 
ambient surroundings as well as radiant exchange with the heat flux 
gage face. In Eq. (7), q p is the net radiation heat transfer out of the 

graphite plate surface. It consists of two parts, the net radiation heat 
flux from graphite plate to the gage q pK and the net radiation heat 

flux from graphite plate to the air q pa . These heat fluxes are obtained 

as 
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where A is the specific resistance. It is also a function of temperature. 
As shown in Fig. 6, it decreases from the room temperature at 25 C 
until around 500° C and thereafter it increases. At the final 
experimental temperature of I840°C, it has almost the same value as 
the one at room temperature. The temperature dependence employed 
in the numerical calculation is interpolated from the curve provided 
by the vendor and shown in Fig. 6. 



Figun 6: ATJ Graphite Electrical Resistivity vs. Temperature 


q pa ={\-F ps )ae p {T*-T*) (9) 


The shape factor F ps is calculated by doing the radiation 

network problem using a 2-D model introduced by Patankar(1989). It 
changes with the geometry of the graphite plate and the distance 
between the graphite plate and the gage. 

In Eq. (4) and Eq. (5), h cx is the average natural convection 
coefficient for the horizontal plate surface 


h rv = h,. r 




exp B h - 1 


(10) 


where the blowing parameter B h for heat transfer is defined as 
n u Cp l +n 2lS Cp 2 

B h= — 0 (II) 

Kx 

In Eq. (10) heat transfer coefficient for zero mass transfer, h cx ° , 
for the heated plate facing up is obtained as 


The boundary conditions employed in the model are: 

y = H p ; 

dT u 

- k p — = h cx u (T - T e ) + a E P (T 4 - r * ) ^ - m x U A H c (4) 

y = 0 ; 

dT (I 
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x^0 ; 
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and for the heated plate facing down as 

0.82 (Gr Pr) 1 ' 5 
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where the Grashof number is 

Gr = p AT g l? /v 2 

where L is the length of the surface. 
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In Eq. (6) and (7), h cy is the average natural convection heat 
transfer coefficient for the vertical plate surface and is given by 
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The heat transfer coefficient for zero mass transfer is 


(14) 


0 0.68 + 0.67- (Gr Pr (p) 


1/4 


777, 


(15) 


For the experimental conditions, at steady state ( t~1840° C ), Gr 
is less than 10 9 . Laminar boundary layer occurs on the graphite plate 
surface. The Prandtl number function <p is defined as 


<? = [! + (- 
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9/16 j-16/9 
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where D ]2 is the binary diffusion coefficient. 


3.2 Numerical model for the heat flux gage 

The heat flux gage is represented by a solid copper block of 
diameter (i.e., height in the 2-D model), D, and thickness, TH^ as 
shown in figure 5. Boundary conditions on the front face and upper 
and lower surfaces are free convection of air and radiation. Forced 
convection of water is the boundary condition on the back face. The 
two dimensional governing energy equation for the heat flux gage is 


where k s is the thermal conductivity of the gage. 
The boundary conditions are: 


(23) 


As a result of combustion, 3.28x1 0 7 J/kg of heat is released. The 
mass transfer rates m are 
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where _ — i s the mass transfer blowing factor and p [fn is the 

Plm 

mass driving force. It is obtained as 
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where m l v is the species 1 mass concentration on the surface of 

the graphite plate, denoting 1 as oxygen. Since the chemical kinetics 
are so rapid, the reaction is diffusion controlled. The gas mixture at 
the graphite plate surface is in the equilibrium with the solid carbon. 
The equilibrium data indicates that the concentration the oxygen at the 
surface is essentially zero. That is =0. Where m i e is the species 

concentration in the ambient air. It equals to 0.231 . The ratio between 
oxygen mass flux transferred to the graphite plate surface and the total 

mass transfer rate is fl I v / in . 

By using the analogy between heat and mass transfer, the zero- 
mass-transfer-limit mass transfer conductance are 
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The convective heat transfer coefficients h cxg and h cyf , are 

given by substituting D for H p in Eq. (15) and substituting TH s for 
TH p in Eq. (12) and Eq. (13), respectively. 

In a manner similar to q f) , q g is the net radiation heat transfer 

out of the gage surface. It is also obtained by including two parts 
which are 
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where q np and q ^ are the net radiation heat flux as from heat flux 

gage to the graphite plate and ambient surroundings, respectively. 

Forced convection cooling occurs due to water flow over the 
back face of the gage substrate. Water comes in the tube 1 on one side 
of the gage, circles around the channel inside the gage and then leaves 
from the tube on the other side of the gage. The inlet water 
temperature is 25° C. The water velocity V is assumed to be - Im/s. 
Reynold number corresponding to this velocity is approximately 
9xl0 4 . Thus, for turbulent flow the heat transfer coefficient, h water , 
is obtained as 
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sensitivity, and possible reasons for errors between the experimental 
(30) data and numerical calculations. 


Where DI is the inside diameter of the channel as shown in Fig. 

5. 

3.3 Erosion Calculation 

During the heating of the graphite plate, there is chemical 
reaction which involves burning carbon in the air. The reaction is 
C+0 2 =C0 2 . Since carbon is continuously consumed at the plate 
surface, the thickness and the width of the graphite plate decreases 
with time. A simplification is made in the model by assuming that 
each point on the plate surface undergoes the same average erosion. 
As a result, the shape of the graphite plate for proposes of 
computation always remains square whereas the actual plate assumes 
a tapered shape. 

The mass balance on the graphite plate surface is 

^r(V p p ) = -m ■ A (31) 

at 

where p p is the density of the graphite plate. The volume of the 
graphite plate, V, is equal to L p • H p TH p and the surface area A of 
one side of the plate and is equal to L p H p . Thus the rate of decrease 
of the thickness is obtained as 

-^(77/ ) = — (32) 

dt p p t> 


3.4 Calculation Procedure 

Due to erosion, the shape of the graphite plate changed with 
time. Yet, based on the experiment, the plate temperature doesn’t 
change rapidly, such that at each erosion time step (one second), the 
temperature is quasi steady. In carrying out the numerical 
calculations, at each erosion time step, the following steps were used; 

1. Give the initial temperature of the graphite plate and the heat 
flux gage. 

2. Apply boundary conditions to every equation and obtain new 
temperature. 

3. Use the graphite plate new temperature in the radiation 
boundary condition of heat flux gage, update temperature of the heat 
flux gage and use it for the radiation boundary condition of graphite 
plate. 

4. Repeat step 2 onwards until the results converge to a unique 
value. 

5. From the erosion rate of the graphite plate, calculate the shape 
change due to combustion and thereafter obtain the new shape factor 
between the graphite plate and the heat flux gage to be used for the 
next time step. 

6. Repeat from step 2 to step 5 until the final time. 


4.1 Temperature distribution 

The temperature distribution in the graphite plate and the heat 
flux gage, as computed from the numerical model, are shown in Fig. 

7. The results shown in Fig. 7 were computed at test time = 60 
seconds, the initial time in the analysis. The measured surface 
temperature at the vertical center of the graphite plate was 2066K. 
This temperature is about 4% higher than that predicted from the 
numerical calculation (1986K). 



Graphite Plate 



0 . 00 . 10 . 20.3 


Thickness (cm) 
Heat Flux Gage 


Figure 7: Graphite Plate and Heat Flux Gage Temperature 
Distribution 


4. RESULTS AND DISCUSSIONS 

Comparison is made between the experi menial data and the 
numerical calculations. Discussion of the results will include 
temperature distribution, heat flux, graphite plate erosion, grid 
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Figure 8: Comparison of Experimental and Numerical Calculated 
Gage Heat Flux 

4.2 Heat flux 

Given a nominal experimental condition for heat flux of 477 
kW/m\ applied for lOmins, a comparison is made between the 
measured heat flux vs time and the numerical lical (lux vs time. Due 
to the combustion, the graphite plate become thinner and thinner. The 
net radiation heat flux from the graphite plate reaching the gage 
decreases with time as shown in Fig. 8. The difference between the 
experimental and numerical prediction is less than 1 %. 


4.4 Grid sensitivity 

The numerical model described in section 3 was first tested for 
its accuracy by comparing the calculated gage heat flux with the 
experimental data using the 5x5 mesh size. The error was -5%. For 
testing the grid resolution, the computations were made by using the 
mesh densities from 7x7 to 107x107. The numerical results are 
plotted in Fig. 10. From Fig. 10, it can be seen that as the number of 
mesh points increase, the heat flux decreases and is approaching an 
asymptotic value for large number of mesh points. As the difference 
between the 37x37 and the asymptotic heat flux values is 
approximately 1%, the computations in this study were performed 
with 37x37 mesh to save the computation time with sufficient 
accuracy of numerical result. 



Grid Poinls Per Axis (square grid) 


4.3 Graphite plate erosion 

Since the combustion of the carbon occurs at the surface of 
the graphite plate, carbon is burned into C0 2 , and the plate become 
thinner with time. It was found that in 10 minutes, the height of the 
plate decreased from 4.14 cm to the final measurement 3.81 cm 
whereas the initial thickness decreased from 0.406 cm to 0.282 cm. 
Both of the experimental and the numerical results are shown in Fig. 
9. 
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Figure 9: Change of the Graphite Plate Height and Thickness 



Figure 10: Grid Sensitivity 

4.5 Possible reasons for the discrepancies found between 
experimental results and numerical computations 

From the comparisons discussed in section 4.1 , 4.2, 4.3, several 
reasons for the difference between the numerical prediction and the 
experimental data are possible. These include 
1 . Experiment error 

The various measurements taken during the course of the 
experiment have uncertainties associated with them. Sources of these 
uncertainties include the readings of the voltmeter, the current meter, 
and the voltage transducers, as well as the electric current 
measurement coil and the optical pyrometer. 

The uncertainty estimates for various measurements are 
summarized in Table 1. The total uncertainty in a measurement is the 
square root of the sum of the squares of the component uncertainties. 
Uncertainties in the heat flux gage calibration are not included since 
the calibration system characterization project is attempting to define 
those uncertainties. 

Uncertainties in the plate dimension measurements are as 
follows: ±0.25 mm for thickness and 2.54 mm for height. A small, 
but unknown, amount of erosion occurs prior to the beginning of the 
analysis. This erosion occurs during a short heating of the plate to 
verify proper clamping of the plate and during the transient startup of 
each experimental run. 
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Table 1: Experimental Measurement Uncertainty 



Note: All source uncertainties are derived 
from manufacturers data. 


2. The effect of parameters 

The properties used in the numerical model came from a variety 
of sources. Material properties of the ATJ graphite were obtained 
from the material manufacturer. Properties for the copper block 
representation of the heat flux gage were taken from a handbook. 
Emissivity of the graphite plate and black paint which coated the heat 
flux gage face were obtained using an emissivity measuring optical 
pyrometer described by Matthews(1989). Most properties were 
obtained as functions of temperature. However, some properties, such 
as the emissivity of the paint on the gage face, could only be obtained 
for room temperature. 

Errors are introduced into the numerical calculation since the 
manufacturer and handbook data represent typical values as opposed 
to measurements for the actual pieces of material used. Properties 
known only at room temperature are likely to vary somewhat with 
temperature thereby introducing another error into the calculation. 

Electrical current passing through the plate is the primary energy 
source in the model and was obtained from measurements. These 
measurements have an uncertainty interval associated with them. 
This uncertainty is another source of error in the numerical heat flux 
computation. Table 2 summarizes the effect variations in several 
parameters have on the numerical results. 


Table 2: Heat flux changes caused the uncertainties of various 
parameters 


Emissivity 

Heat 

Flux 

Current 

Heat 

Flux 

Thermal 

Conductivity 

Heat 

Flux 

±0.1 

11 % 

±24 Amp 

4% 

±5 

.04% 


anymore and this can affect the shape factor then finally affect the 
estimation of the heat flux received by the gage. 

The numerical analysis currently models the heat flux gage as a 
solid piece of copper. In reality, the circular foil heat flux gage used 
in the experiment has a far more complex structure. This 
simplification does not effect incident heat flux computations but will 
impact absorbed heat flux calculations. The heat flux data presented 
in this paper are incident heat flux. The detail of the actual heat flux 
gage will be incorporated into future versions of the numerical 
analysis. 


5. CONCLUSIONS 

A two-dimensional finite difference model of a flat plate heat 
flux gage calibration system has been developed. The numerical 
solution provides temperature distributions through the thickness of 
the flat plate heating element and the heat flux gage, the heat flux at 
the gage face, and the shape change of the flat plate with time. Initial 
grid refinement has been performed and has reduced grid error to less 
than 0.05% between successive grid sizes. Preliminary comparison to 
experimental data has also been performed. Differences between 
numerical and experimental data are less than 7%. 

Additional work needs to be performed in order to identify and 
quantify the effect of various error sources present in the calibration 
system. This work will incorporation of additional detail of the heat 
flux gage structure, and extension of the model to 3 dimensions. The 
numerical analysis will eventually be extended to uncooled heat flux 
gages operating at high (850° C) temperature. 
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3. The limitation of the numerical model 

In reality, the graphite plate shape varies in all the 3-dimensions. 
The incident heat flux can be estimated more accurately only when 
three-dimensional computation is performed. Also, even for the 2-D 
model, the graphite plate burning rate would be different at different 
location on the surface since there is low mass transfer rate at the top 
of the plate, and high mass transfer rate at the bottom of the plate. 
Thus, the final shape of the cross section does not remain rectangular 
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